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A  COMPUTATIONAL  MODEL  FOR  CONTACT  STRESS  PROBLEMS 
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Virginia  Polytechnic  Institute  and  State  University 
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SUMMARY 

A  mixed  finite  element  model  based  on  the  incremental  nonlinear 
theory  of  solid  bodies  is  developed  and  its  numerical  accuracy  is 
evaluated  by  applying  the  model  to  a  number  of  contact  problems.  The 
theoretical  basis  of  the  finite  element  model  is  formulated  using  the 
updated  Lagranqian  formulation.  Both  displacement  and  mixed  models  are 
described,  but  emphasis  is  placed  on  the  development  and  application  of 
the  mixed  model,  which  contains  the  displacements  and  stresses  as  the 
nodal  degrees  of  freedom. 


INTRODUCTION 

In  the  linear  description  of  the  motion  of  solid  bodies  it  is 
assumed  that  the  displacements  are  infinitely  small  and  that  the 
material  is  linearly  elastic.  In  addition,  it  is  also  assumed  that  the 
nature  of  the  boundary  conditions  remains  unchanged  during  the  entire 
deformation  process.  These  assumptions  imply  that  the  displacement 
vector  u  is  a  linear  function  of  the  applied  load  vector  F  ,  i.e.,  if 
the  applied  load  vector  is  a  scalar  multiple  oF  then  the  corresponding 
displacements  are  au  . 

The  nonlinearity  in  solid  mechanics  arises  from  two  distinct 
sources.  One  is  due  to  the  kinematics  of  deformation  of  the  body  and 
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the  other  from  the  constitutive  behavior  (i.e.,  stress-strain 
relations).  The  analyses  in  which  the  first  type  of  nonlinearity  is 
considered  are  called  geometrical ly  nonlinear  analyses,  and  those  in 
which  the  second  type  are  considered  are  called  materially  nonlinear 
analyses.  The  geometrically  nonlinear  analysis  can  be  subclassified 
according  to  the  type  of  nonlinearities  considered.  Two  such  cases 
are:  (i)  large  displacements,  large  rotations,  but  small  strains,  and 
(ii)  large  displacements,  large  rotations  and  large  strains.  In  the 
first  case  it  is  assumed  that  rotations  of  line  elements  are  large,  but 
their  extensions  and  changes  of  angles  between  two  line  elements  are 
small.  In  the  second  case  the  extension  of  a  line  element  and  angle 
changes  between  two  line  elements  are  large,  and  displacements  and 
rotations  of  a  line  element  are  also  large. 

The  objective  of  the  present  study  is  two-fold:  First,  to  describe 
the  equations  that  govern  the  nonlinear  behavior  of  solids,  and  second, 
to  develop  finite  element  models  for  the  analysis  of  the  nonlinear 
behavior  of  solid  bodies  in  contact. 

INCREMENTAL  EQUATIONS  OF  MOTION 

Consider  the  motion  of  a  body  in  a  fixed  Cartesian  coordinate 
system.  Suppose  that  the  body  can  experience  large  displacements,  large 
strains  and  nonlinear  mechanical  (i.e.  constitutive)  behavior.  We  wish 
to  determine  the  configuration  of  the  body  for  different  times  and 
loads.  The  formulation  to  be  described  assumes  that  the  configurations 
of  the  body,  c^,  C2  .  .  .  ,  due  to  loads  Pj,  ,  P^, 

respectively  are  known  and  seeks  the  configuration  at  time  Thus, 

in  the  present  formulation  we  follow  the  body  as  it  deforms  from  the 
initial  configuration  to  the  final  configuration.  This  type  of 


description  is  called  the  Laqranqian  (or  material)  description,  which 
differs  from  the  Eulerian  description  used  in  fluid  mechanics 
problems.  In  the  Eulerian  description  of  motion,  instead  of  following 
the  body  (or  a  fixed  collection  of  particles  constituting  the  body),  the 
motion  of  particles  passing  through  a  fixed  control  volume  is  determined 
for  various  times.  The  Lagrangian  description  is  a  natural  one  for  a 
solid  body  because  one  is  more  interested  in  the  deformation  of  the  body 
than  in  the  changes  that  are  taking  place  in  the  control  volume  that  was 
occupied  initially  by  the  body.  On  the  other  hand,  the  Eulerian 
description  is  a  natural  one  for  fluid  flow  problems,  because  there  one 
is  interested  in  the  velocity,  pressure  and  stress  fields  in  a  fixed 
control  volume  without  focusing  attention  on  fluid  particles  that  enter 
and  leave  the  control  volume. 

In  the  Lagrangian  description  of  motion  all  variables  are  referred 
to  a  reference  configuration,  which  can  be  the  initial  configuration  or 
any  other. convenient  configuration.  The  description  in  which  all 
variables  are  referred  to  the  current  configuration  is  called  the 
updated  Laqranqian  description  [1-4]  and  the  one  in  which  all  variables 
are  referred  to  the  initial  configuration  is  called  the  total  Laqranqian 
description  [5-8]. 

The  equations  of  the  Lagrangian  incremental  description  of  motion 
can  be  derived  from  the  principles  of  virtual  work  (i.e.,  virtual 
displacements,  virtual  forces  or  mixed  virtual  displacements  and  forces) 
[9-13].  Since  our  ultimate  objective  is  to  develop  the  finite  element 
models  of  the  equations  governing  a  body,  we  will  not  actually  derive 
the  differential  equations  of  motion  but  utilize  the  virtual  work 
statements  to  develop  the  finite  element  models. 


Displacement  Formulation 

The  displacement  finite  element  model  is  based  on  the  principle  of 
virtual  displacements  (see  Reddy  [14]).  The  principle  requires  that  the 
sum  of  the  external  virtual  work  done  on  a  body  and  the  internal  virtual 
work  stored  in  the  body  should  be  equal  to  zero: 


!  dv  -  m'f)  =  0 


(1) 


where 


^ij 


2®ij 


•f . 

1 


■t.  = 


6(^F)  =  f  ^f.  6U.  dV  +  J  ^t.  6U.  dS 

the  Cartesian  components  of  the  Cauchy  stress  tensor  in 
configuration  C2  at  time  (t  +  At) 
occupying  the  volume 

the  Cartesian  components  of  the  infinitesimal  strain 
tensor  associated  with  the  displacements  u^-  in  going  from 
configuration  c-^  at  time  t  to  C2  at  time  (t  +  At)  : 

3U.  3U. 


p  =1  ( _ I  +  . 

2^ij  2  '‘3X  .  3x.^ 

J  ’ 


(2) 


Cartesian  components  of  a  point  in  configuration  C2. 
Cartesian  components  of  the  body  force  vector  measured 
in  C2 

Cartesian  components  of  the  surface  stress  vector 
measured  in  C2 


Here  5  denotes  the  variational  symbol  (i.e.,  6u.  denotes  the  virtual 
displacement  in  u^)  and  dV  and  dS  denote  the  volume  and  surface  elements 
in  the  configuration  over  which  the  integrals  are  defined. 
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Equation  (1)  is  not  so  useful  in  its  present  form  because  the 
integrals  are  defined  over  the  volume  V2  and  surface  $2  of  the 
configuration  C2,  which  is  yet  unknown.  In  the  linear  analysis,  it  is 
assumed  that  the  configuration  of  the  body  remains  unchanged  and 
therefore  Eq.  (1)  applies  to  the  initial  (undeformed)  configuration. 

The  fact  that  the  configuration  of  the  body  changes  continuously  in  a 
nonlinear  analysis  requires  us  to  use  appropriate  measures  of  stress  and 
strain  and  their  interrelationship  (i.e.,  constitutive  equations)  so 
that  Eq.  (1)  can  be  used  to  calculate  the  configuration  c^.  In  order 
to  compute  the  current  configuration  C2  (often,  the  displacements, 
strains  and  stresses  in  C2)  from  the  knowledge  of  applied  forces  and 
displacements  and  known  previous  configurations  Cj^,  we  must  make  some 
assumptions.  A  description  of  the  procedure  based  on  the  updated 
Lagrangian  approach  is  given  below. 

The  coordinates  of  a  general  point  in  and  Cj  and  C2  are  denoted 
by  (X°,  X^,  X°),  (X^,  X^,  X^),  and  (x^,  x^,  x^),  respectively.  The 
displacements  of  a  general  point  in  are  denoted  by  (^u^,  ^u^,  ^u^). 

In  C2  they  are  given  by 

^u.  =  ^u.  +  u.  ,  i  =  1,2,3  (3) 

where  u^  are  the  components  of  the  displacement  vector  from  c-^  to  C2. 

During  the  motion  of  the  body,  its  volume,  surface  area,  density, 
stresses  and  strains  change  continuously.  The  stress  measure  that  we 
shall  use  is  the  second  Piola-Kirchhoff  stress  tensor.  The  components 
of  the  second  Piola-Kirchhoff  stress  tensor  in  will  be  denoted  by 
S^-j.  To  see  the  meaning  of  the  second  Piola-Kirchhoff  stress  tensor, 
consider  the  force  dF  on  surface  dS  in  C2-  The  Cauchy  stress 


tensor  t  is  defined  by 


(n  •  t)  dS  =  dF  (4a) 

where  n  is  the  unit  normal  to  dS  in  02-  Note  that  the  Cauchy  stress  is 
the  force  per  deformed  area  (i.e.  measured  in  C2)  and  referred  to  c^- 
The  second  Piola-Kirchhoff  stress  tensor  at  time  t  +  At  referred  to  Cj 
is  defined  by 

(%  •  l5>  ‘‘^0  = 

where  n  denotes  the  unit  normal  to  the  surface  element  dSg  in  Cj^.  The 
force  dF^  is  related  to  dF  by 

J’^  •  dF  (4c) 

where 

1  3X  T 

r'  ■  (if) 

From  the  definition  it  is  clear  that  the  second  Piola-Kirchhoff  stress 

is  measured  in  C2  but  referred  to  Cj^.  It  can  be  shown  that  (see 

2  2 

Malvern  [15])  that  the  components  S. .  and  x. .  are  related  according 

*  V  ’  J 

to 


p  3X.  „  3X  . 

2s  =  _o  _ 1.2^  , 

1  ij  p  3Xj^  mn  3X^ 


(5a) 


2 

T 


ij 


3X. 

1 


(5b) 


where  p^  denotes  the  density  in  Cj^  and  p  is  the  density  in  C2.  The 

second  Piola-Kirchhoff  stress  tensor  is  symmetric  whenever  the  Cauchy 

2  2  2 

stress  tensor  is  symmetric.  Note  that  -S..  = 


Similarly,  the  Green- Lagrange  strain  tensor  and  the 
infinitesimal  strain  tensor  e^j  are  related  by 

tp  ^  m n  / c\ 

Tij  '  aX.  aXj  2®mn 

It  is  also  important  to  note  that  the  2nd  Piola-Kirchhoff  stress 
tensor  is  energetically  conjugate  to  the  Green-Lagrange  strain  tensor 
and  the  Cauchy  stress  is  energetically  conjugate  to  the  infinitesimal 
strain  tensor.  In  other  words,  we  have 


Substituting  the  eguality  (7)  into  Eg.  (1),  we  obtain 

“  -  !  Kj  («) 

''l 

Next  we  use  the  incremental  decompositions  of  the  stresses  and  strains: 

2c  -  1  ^  c 

l^ij  ■  ^-j  l^ij 


iE.,  =  iS-.  + 

1  ij  1  ij  1 


where 


jS^j  =  incremental  components  of  2nd  Piola-Kirchhoff  stress  tensor 

le.- ,•  =  (incremental)  components  of  the  infinitesimal  strain  tensor 
,  au.  3u. 

2  ^X.  aX.^’ 

1  au„ 

i  _ m  _ m  n n^ 

I'^ij  -  2ixTix- 

Recall  that  u^-  is  the  i-th  displacement  component  of  a  generic  point  in 
(in  going  from  c-^  to  c^).  Substituting  Eg.  (9)  into  Eg.  (8),  we 
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U 

or 

!,  ihi  «(i^j  *  I'-ij)  <'''  *  I  ‘'ij«(i'ij)<iv 

■  -  J'„  ‘<i'u>‘‘''  " 

Linearize  the  equations  by  assuming  that 


l^ij  l^ijrs  ®rs’  ^  l^ij  '^l^ij 
We  obtain  the  approximate  equation  of  equilibrium, 


!  I'^ijrs  l^rs  O''  !  ‘'U  «(l''ij)1'' 

Vi 

-  "A,  Sj  (‘3) 


This  linearization  can  be  interpreted  as  a  representation  of  the 
nonlinear  curve  between  two  consecutive  load  steps  by  a  linear  line 
segment. 


Mixed  Formulation 

A  mixed  (or  stationary)  virtual  work  statement  that  treats  the 
displacements  and  stresses  as  independent  (dependent)  variables  can  be 
derived  from  Eq.  (13)  as  follows.  First,  we  note  that  Eq.  (13)  is  the 
first  variation  of  the  functional 

"  ’  J  I'^ijrs  iV  l"ij  "  J'v  ‘’Ij'l'ij  *  04) 

'^1  ''l 

Clearly,  .1  denotes  the  (approximate)  total  potential  energy  of  the  body 
in  configuration  C2  but  referred  to  configuration  c-^,  and  Eq.  (13)  is  a 


statement  of  the  principle  of  the  minimum  total  potential  energy.  In 

Eqs.  (13)  and  (14)  it  is  understood  that  .e..  and  ,n. •  are  defined  in 

»  •  3  '  3 

terms  of  u^-  [see  Eq.  (10)]. 

In  the  mixed  formulation,  we  begin  with  the  expression  for  the 
approximate  total  potential  energy  n  and  introduce  the  stresses  as 
variables  using  the  Lagrange  multiplier  method.  We  treat  the  strain- 
displacement  relations 


as  constraints  and  let  the  increment  in  the  2nd  Piola-Kirchhoff  stress 
tensor  act  as  the  Lagrange  multiplier.  Including  this  in  the 
variational  form  results  in  the  modified  functional  iT|^,  which  is  given 
by 


l^ij4®ij 


(16) 


where  n  is  given  by  Eq.  (14).  Here  we  note  that  the  sign  of  the 
Lagrange  multiplier  term  is  generally  arbitrary,  and  is  selected  as 
negative  to  obtain  the  correct  form  of  the  equations  that  follow. 

The  variational  principle  (often  called  the  Hel 1 inger-Rei ssner 
variational  principle)  is  given  by  setting  5tT|^  =  0.  This  gives 


‘’R  *  ^  l^J  ‘’ij®  I'lj  *  l®k.  -  l^ij  l=ij) 

''l 

iS,j.u,_jldV  -  «(2f)  .0  (17: 


We  next  write  the  expressions  for  the  strain  energy  density  and  the 

*■ 

complementary  strain  energy  density  due  to  the  incremental 
displacements  as 


displacements  as 


and 


^0  "  2  ^ijki  l®ij 


*  1 

U  =  4  D.  ,S.  .  ,S, 

0  2  ijki  1  ij  1  ka 


(18) 


where  are  the  components  of  the  compliance  tensor, 

write 


We  may  then 


U  =  U  -  .S. .  ,e. . 
0  0  1  ij  1  1J 


(19) 


and,  taking  the  variation  of  this  expression  gives 


*  1 

5U  =  6(4  C.  ,e.  .  .e.  -  ,S.  .  ,6.  .) 

0  '2  ijkn  1  ij  1  ka  1  Tj  1 


(20) 


If  we  write  the  complementary  strain  energy  density  in  terms  of  the 

2nd  Piola-Kirchhoff  stress  components  as  we  did  in  Eq.  (18)  we  may  write 

* 


-  6U^  =  - 


3U 

*  l^ij  "  "  l^ka^  l^ij 


(21) 


Equating  the  expressions  in  Eqs.  (20)  and  (21)  and  substituting  into  Eq. 
(17)  yields  the  two  approximate  equilibrium  equations 


;  u,_j.(iS,j)dv- j  .)dv.  0 


(22) 


These  two  equations  are  analogous  to  Eq.  (13)  for  the  displacement 
formulation. 

In  the  next  section  we  discuss  the  finite  element  models  of  Eqs. 
(13)  and  (22). 
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2^ 


FINITE  ELEMENT  MODELS 


Displacement  Model 


Here  we  construct  the  finite  element  model  of  Eq.  (13)  for  the  two- 
dimensional  case.  We  assume  that  the  displacement  components  are 
interpolated  by  expressions  of  the  form  (see  Reddy  [16]) 

n  . 

u^-(Xj,X2)  -  u.  ijj  (Xj,  X2)  (23) 

J  ^ 

where  o'?  denote  the  value  of  u^  at  the  j-th  node,  and  denote  the 
interpolation  functions.  Substituting  Eq.  (22)  into  Eq.  (13)  we  obtain 

([k'-]  +  (K°l)  {^}  =  {F^}  -  {F°}  (24) 


where 


[K*-]  =  t  f  [8^]^  [C]  [b'-]  dA  ,  t  =  thickness 
A, 


^11  ^12  ^16 
[C]  =  Cj2  ^22  ^26 

^16  ^26  ^66 


n,i  °  '^2,1  ° 

(B^l  =0  ip.  «  0  ipp  « 

(3  X  2n)  , 

n,2  n,l  '^2,2  '^2,1 


[K°/  =  t  /  (8'^]'  (rl  IB^]  dA 

A, 


ip_  1  0 

n,l 

0  )p_  ~ 

n,Z 

n  y  C  9  i 


IB^l 

(4  X  2n) 


11 

^12 

u 

u 

12 

^22 

0 

0 

0 

0 

^11 

M2 

0 

0 

^12 

^22 

{F*-}  =  t  /  M'  {f}dA.  {F°}  =  t  J  [B*-)'  {T}dA 


f'^'i  0  0  .  .  .  on  (m)  (^1  ) 

(2  <  2n)  Lo  *10  *2  ...  0  If  i  |  j 


It  is  understood  that  is  computed  using  the  Almansi  strain 

'  xj 

tensor, 

^ij  ”  l^ijkm  l^km 


,  3Ut  3U„  3U„  3U„ 

—  A  /  4.  ni  in  m\ 

l^km  ■  2  ^3X„  3X.  ■  3X.  3Xj 

m  k  km 


where  the  displacements  are  now  the  total  displacement  components 
measured  in  the  current  configuration.  Also,  since  Eq.  (13)  is  a 
linearized  version  of  Eq.  (11),  the  error  introduced  into  the 
calculation  of  the  displacements  u.j  between  configurations  can  drift  the 
solution  away  from  the  true  solution  (especially  if  the  load  steps  are 
large).  Therefore,  a  correction  should  be  made  to  the  displacements  at 
each  load  step.  This  can  be  done  as  follows:  The  solution  {a}  of  Eq.  (24) 
allows  us  [with  the  aid  of  Eq.  (23)]  to  compute  the  total  displacements 
at  time  (t  +  at)  , 

\  =  ‘>i,  *  u, 


which  can  be  used  to  compute  the  strains  and  stresses  (with  appropriate 
constitutive  equation)  at  time  t  +  at.  By  the  principle  of  virtual 
displacements,  the  true  displacements,  strains  and  stresses  at  any  time, 
say  at  time  t  +  at,  are  such  that  the  internal  virtual  work  is  equal  to 
the  external  virtual  work  done.  Since  u^-  (hence  the  strains  and 
stresses  computed  from  them)  are  approximations,  there  will  be  an 
imbalance  between  the  internal  and  external  virtual  works  performed  on 
the  body.  This  imbalance  can  be  minimized  by  updating  the  internal 
virtual  work  through  an  iteration  (for  a  fixed  system  of  loads  and 
time);  the  iteration  is  continued  until  the  imbalance  is  reduced  to  a 
preassigned  value  (i.e.,  a  convergence  limit).  For  example,  the 
displacement  increment  at  the  (r  +  l)st  iteration  is  calculated  from  the 
equations 

({k'-I  +  (K^^l)^  ^  1  =  {fS  -  {F"}^  (27) 

wherein  [t]  and  {t}  are  calculated  using  the  displacements, 

"  (“I'r 

Equations  (27)  and  (28)  correspond  to  the  Newton-Raphson 
iteration.  If  the  left  hand  side  is  not  updated  during  the  iteration, 
the  iterative  scheme  is  known  as  the  modified  Newton-Raphson 
iteration. 

Mixed  Model 

We  next  construct  the  finite  element  model  of  Eqs.  (22)  for  the 
two-dimensional  case.  We  begin  by  assumming  independent  approximations 
of  the  displacements  and  stresses  of  the  form 


(29) 


u.(Xj,X2) 


z  u}«n.(x,,x„) 
j=l  1  J  ^ 


l^ij  l^ij'*'k^^l’^2^ 

Substituting  these  expressions  into  the  two  equilibrium  equations  in  Eq 
(22)  gives 

(30) 

where 

(K^^l  =  tj  [B°]'^lTl[B°ldA 
(K^^l  =  tJ  {^°]''’[Dl[f“]dA 

(K^2)  =  tJ  fB‘-]'^h°]dA  (31) 

(f'-}  =  tf  (fl'^{f}dA 
=  tJ  (BLl^lTldA 


All  of  the  other  necessary  matrices  have  been  defined  previously  during 
the  development  of  the  displacement  finite  element  model. 

As  with  the  displacement  model,  Eq.  (30)  is  solved  repeatedly  until 
the  force  imbalance  is  reduced  to  a  fixed  tolerance.  This  amounts  to 
measuring  the  percentage  of  the  displacement  and  stress  increments 
relative  to  the  total  solution  vector.  When  the  increase  in  the 
displacements  and  stresses  has  been  reduced  to  below  a  very  small 
percentage  of  the  total  solution,  the  state  of  equilibrium  has  been 
obtained  for  the  given  load  step,  and  the  load  may  be  increased  or  the 
analysis  may  be  terminated. 

The  solution  of  Eq.  (30)  allows  us  to  compute  the  total 
displacements  and  Cauchy  stresses  at  time  (t  +  At).  To  compute  the 
displacements  we  simply  use  the  equation 

^u .  =  ^u .  +  u . 

1  ^  1 

In  the  mixed  formulation,  there  is  no  need  to  compute  the  Cauchy 
stresses  using  the  Almansi  strains  as  was  required  using  the 
displacement  formulation.  Since  the  increments  of  the  2nd  Piola- 


Kirchhoff  stress  tensor  are  computed  as  nodal  variables,  we  simply  use 
our  incremental  decomposition  of  the  stress  given  in  Eq.  (9)  as 

iSlj  -  T.^.  +  jS.j 

When  the  increment  in  the  ,S. .  terms  is  reduced  to  be  within  the 

L  1  J 

2  2 

required  tolerance,  we  have  tS..  =  t... 
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PLANAR  ELASTIC  CONTACT  PROBLEMS 

In  this  section  we  discuss  several  techniques  for  the  analysis  of 
two-dimensional  elastic  contact  problems.  Such  problems  have  a  host  of 
computational  difficulties  since  the  region  of  contact  is  typically  not 
known  in  advance,  nor  are  the  regions  of  relative  stick  and  slip  between 


the  two  contacting  bodies  due  to  the  presence  of  friction.  Most  current 
algorithms  that  solve  contact  problems  are  relatively  complex,  and  use  a 
number  of  iterative  schemes  to  account  for  the  changing  boundary 
conditions  and  regions  of  contact. 

The  nonlinear  mixed  finite  element  formulation  forms  the 
cornerstone  of  the  methods  described  below.  The  displacement  finite 
element  method  has  been  used  almost  exclusively  in  previous  numerical 
analyses  of  contact  problems.  Mixed  elements  provide  the  immediate 
advantage  of  computation  of  stresses  as  nodal  variables,  which  is  ideal 
for  contact  problems  since  the  stresses  may  be  obtained  precisely  on  the 
contact  boundary. 

Two  different  and  independent  methods  are  described  for  the 
analysis  of  the  specific  problem  of  a  pin-loaded  plate.  The  methods  are 
described  below  and  several  examples  are  presented  for  the  first  of  the 
two  techniques. 


Rigid  Pin  Contact  Algorithm 


In  general,  contact  problems  involve  two  or  more  elastic  bodies 
pressing  against  one  another  under  external  load.  In  the  case  of  a  pin- 
loaded  plate,  these  two  bodies  are  the  plate  and  the  pin.  If  the 


assumption  of  a  rigid  pin  is  used,  the  analysis  is  simplified 
considerably.  This  assumption  eliminates  the  need  to  analyze  the  pin, 
which  not  only  provides  a  known  point  of  reference  for  the  ensuing 
contact  (i.e.  the  surface  of  the  pin),  but  also  drastically  reduces  the 
resulting  global  finite  element  system  of  equations  since  there  is  no 
need  to  discretize  the  domain  of  the  pin. 


The  assumption  of  a  rigid  pin  is  reasonable  if  the  modulus  of 
elasticity  of  the  pin  is  much  higher  than  that  of  the  plate.  Analytical 


studies  have  also  shown  that,  in  the  contact  analysis  of  composite 
plates,  pin  elasticity  is  not  an  important  variable  and  has  a  relatively 
small  effect  on  the  resulting  stress  distributions  [17]. 

One  simple  and  effective  method  for  analyzing  thin,  orthotropic, 
pin-loaded  plates  was  developed  originally  by  Wilkenson  (18)  and  was 
later  refined  by  Rahman  (19)  to  capitalize  on  the  computational 
advantages  that  arise  from  the  rigid  pin  assumption.  The  method  uses 
three  separate  iteration  steps  to  account  for  the  incremental  load 
level,  the  contact  process,  and  the  effects  of  friction.  Both  the 
original  and  refined  schemes  used  displacement  finite  elements.  In  the 
load  step  iteration,  the  solution  for  a  given  load  increment  was  treated 
as  a  linear  analysis,  i.e.  the  equations  of  linear  elasticity  were 
used.  The  stresses  in  the  plate  were  computed,  for  a  given  load  level, 
using  the  original,  undeformed  configuration  of  the  plate  along  with  the 
final  displacement  vector. 

The  use  of  displacement  elements  in  the  analysis  of  con^'ct 
problems  frequently  results  in  very  large  displacement  gradients  at  the 
region  of  contact.  Since  the  required  stresses  are  computed  at  the 
element  interiors  and  are  then  extrapolated  to  the  contact  boundary, 
some  type  of  stress  smoothing  is  often  necessary  using,  for  example,  a 
local  least  squares  routine  [20)  or  iterative  improvement  on  the 
averaged  nodal  stresses  (21).  Using  mixed  elements,  this  is  not 
necessary  since  the  stresses  are  computed  as  nodal  variables  and  no 
postcomputation  is  necessary  to  modify  the  resulting  nodal  stresses.  It 
is  for  this  reason  that  mixed  elements  would  appear  to  be  advantageous 
over  displacement  elements  for  contact  problems. 


The  refined  algorithm  developed  by  Rahman  [19]  uses  a  mixed  polar- 
Cartesian  coordinate  system  to  fix  the  proper  displacements  of  the  nodes 
of  the  plate  that  have  come  in  contact  with  the  pin.  An  iterative 
scheme  is  used  to  ensure  that,  for  a  given  load  step,  all  nodes  that 
have  come  in  contact  remain  In  contact.  In  other  words,  after  every 
iteration  the  positions  of  all  contact  nodes  are  corrected  in  the  radial 
direction  so  that  they  remain  on  the  surface  of  the  pin.  If  the 
resulting  shearing  stress  for  a  given  contact  node  is  larger  than  the 
induced  radial  stress  multiplied  by  the  nodal  coefficient  of  friction, 
the  node  is  considered  to  be  sliding,  and  it  may  subsequently  move  in 
the  tangential  direction  of  the  pin.  Otherwise,  the  node  is  considered 
to  be  sticking  to  the  pin  due  to  friction,  and  is  fixed  to  an 
interpolated  position  on  the  pin  for  the  remainder  of  the  analysis. 

This  iterative  procedure  is  repeated  until  the  sum  of  the  load  steps  has 
reached  the  required  load  level.  The  details  of  the  this  method  are 
more  completely  described  in  reference  [19]. 

Elastic  Pin  Contact  Algorithm 

The  assumption  of  a  rigid  pin,  which  is  reasonable  for  cases  when 
the  two  bodies  in  contact  have  a  very  large  difference  in  modulus  of 
elasticity,  is  not  usually  valid  for  contact  between  two  generic 
bodies.  The  algorithm  described  in  the  previous  section,  though  useful 
for  certain  problems,  was  mainly  developed  to  demonstrate  the  use  and 
accuracy  of  the  mixed  finite  element  method  for  the  analysis  of  contact 
problems.  For  general  problems  involving  arbitrary  bodies,  it  is 
necessary  to  revise  the  analysis  to  include  the  effects  of  pin 
elasticity,  which  may  be  significant  if  the  two  bodies  are  of  similar 
constitution. 


18 


To  handle  the  complications  arising  from  contact  and  the  presence 
of  friction  between  two  elastic  bodies,  we  add  a  Lagrange  multiplier 
contribution  to  our  original  expression  in  Eq.  (16),  which  will 
represent  the  total  potential  of  the  contact  forces  acting  at  the  nodes 
on  the  contact  boundary.  In  addition,  the  kinematics  of  the  elements  at 
the  contact  interface  must  be  monitored  such  that  the  nodal 
displacements  are  compatible.  We  therefore  invoke  stationarity  of  the 
modified  functional 

k 

’'N  =  "L  - 

where  k  represents  the  number  of  the  contactor  nodes  on  the  boundary  and 
W  represents  the  total  potential  for  a  given  contact  force  acting  at  a 
given  contact  node. 

Taking  the  first  variation  of  the  modified  functional  results  in 
a  system  of  equations  similar  to  Eq.  (30)  with  the  increment  in  the 
contact  forces  added  as  new  nodal  variables  for  the  contactor  nodes. 
These  constraint  equations  may  be  added  to  the  original  stiffness 
matrices  by  means  of  contact  matrices,  which  contain  the  needed 
relationships  between  the  contact  forces  and  the  contact  node 
displacements.  This  idea  was  originated  and  developed  using  a 
displacement  formulation  by  Bathe  and  Chaudhary  [22]. 

The  analysis  proceeds  in  a  manner  similar  to  the  steps  performed 
for  nonlinear  solid  mechanics  problems  in  that  load  increments  are 
applied  to  the  body  or  bodies  and  the  equilibrium  equations  are  solved 
iteratively  until  the  solution  increment  is  within  a  preassigned 
tolerance.  If,  however,  there  is  penetration  of  a  contactor  node  from 
one  of  the  bodies  into  the  domain  of  a  contact  element  from  the  other 


body  during  a  given  iteration,  the  contact  matrices  are  imposed  for  the 
next  iteration  and  the  new  set  of  equations  is  solved  until  the  overlap 
distance  created  during  the  penetration  is  eliminated.  This  results  in 
an  increase  of  the  contact  forces  for  the  nodes  involved  in  the  contact 
as  well  as  a  change  in  the  nodal  displacements  and  stresses  for  these 
nodes.  Inherent  in  the  contact  matrices  are  allowances  for  stick  and 
slip,  and  their  implementation  is  a  function  of  the  contact  stresses  and 
coefficients  of  friction  between  the  two  bodies.  The  coding  of  this 
method  is  underway  and  no  example  problems  have  been  completed  to 
include  in  this  report. 


HYBRID  EXPERIMENTAL/NUMERICAL  TECHNIQUE 

The  primary  difficulty  in  the  analysis  of  contact  problems  is  the 
lack  of  knowledge  regarding  the  region  of  contact  and  the  state  of  stick 
or  slip  between  the  two  bodies.  If  the  regions  of  contact  and  the  state 
of  stick  or  slip  are  known  for  a  given  load  level,  the  analysis  is 
greatly  simplified.  A  hybrid  experimental/numerical  technique  has  been 
proposed  as  part  of  this  research  project  that  combines  the  experimental 
technique  of  moire  interferometry  with  the  numerical  finite  element 
method  to  form  a  hybrid  technique.  The  new  method  uses  the  advantages 
of  each  of  the  methods  to  create  an  accurate  and  powerful  tool  of 
analysis. 

Moire  interfereometry  provides  the  required  sensitivity  that  is 
needed  to  accurately  measure  the  small  in-plane  displacements  of  the 
two-dimensional  pin-loaded  plate  contact  problem.  In  particular,  this 
method  provides  valuable  information  on  the  state  of  the  displacements 
at  the  interface  boundary  between  the  pin  and  the  plate.  This 
eliminates  the  need  for  any  theoretical  assumption  on  exactly  what  is 


occurring  at  the  contact  boundary  that  is  so  common  in  most  numerical 
simulations.  As  mentioned  previously,  it  is  the  complex  state  of 
contact,  stick,  and  slip  that  makes  this  problem  so  difficult  to  model 
in  a  strictly  numerical  method,  and  the  use  of  an  experimental  technique 
greatly  simplifies  the  computational  effort. 

The  exact  displacements  of  a  physical  specimen  may  be  measured  for 
a  sequence  of  increasing  load  steps.  These  displacement  increments, 
along  with  the  loads  applied  to  the  plate,  are  input  as  a  sequence  of 
nonhomogeneous  boundary  conditions  for  the  simulated  problem  analyzed  by 
the  nonlinear  mixed  finite  element  method.  The  resulting  normal  and 
shearing  stresses  around  the  hole  of  the  plate  may  then  be  computed  to 
determine  the  quantitative  behavior  of  the  coefficients  of  friction  as  a 
function  of  the  angular  position  around  the  hole  and  the  load  level. 

The  only  approximations  introduced  are  due  to  the  physical  and 
mathematical  limitations  of  moire  interfereometry  and  the  finite  element 
method.  The  numerically  solved  problem  is  simply  a  special  boundary 
value  problem  with  specified  displacements,  and  no  other  assumptions  are 
involved. 

Since,  in  theory,  the  one  half  of  the  plate  acts  as  the  mirror 
image  of  the  other  half  of  the  plate,  only  one  half  of  the  plate  domain 
is  modeled  for  the  finite  element  analysis.  This  also  drastically 
reduces  the  computational  time  involved.  Preliminary  experimental  data 
have  shown  that  the  displacements  around  the  hole  of  the  plate  are  not 
quite  symmetric,  even  for  an  isotropic  material.  This  is  most  probably 
due  to  the  limitations  of  creating  a  perfectly  symmetric  specimen  and 
applying  a  perfectly  symmetric  load.  The  displacements  are  therefore 


averaged  in  both  of  the  coordinate  directions  of  the  plate  to  yield  one 
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value  for  a  given  point  of  the  contact  boundary  for  the  modeled  half  of 
the  plate. 

The  finite  element  portion  of  the  analysis  currently  awaits  the 
reduction  of  the  experimental  data  for  a  series  of  increasing  load 
steps.  No  examples  of  the  hybrid  procedure  are  included  in  this  report. 

NUMERICAL  EXAMPLES 

In  this  section  we  consider  a  specific  example  to  demonstrate  the 
characteristics  and  differences  of  the  displacement  and  mixed  finite 
element  models  in  nonlinear  analysis. 

Before  considering  the  example  problem,  we  note  an  important  point 
on  the  order  of  approximation  using  the  mixed  elements.  Using  the 
results  of  an  eigenvalue  analysis  of  the  mixed  finite  element  matrices, 
Mirza  and  Olsen  [231  proposed  and  verified  a  completeness  criterion  that 
restricts  the  choice  of  the  order  of  approximation  for  the  displacements 
and  stresses.  The  completeness  criterion  was  given  as: 

The  strains  from  the  stress  approximations 
should  possess  at  least  all  the  strain  modes 
that  are  present  in  the  strains  derived  from  the 
displacement  approx imatons. 

When  this  criterion  is  violated,  the  global  stiffness  matrix  in  the 
mixed  model  will  be  singular  even  after  the  imposition  of  the  boundary 
conditions.  Isoparametric  rectangular  elements  are  used  for  all  of  the 
examples  considered  herein  and,  to  meet  the  requirements  of  the 
completeness  criterion,  only  linear-linear  or  quadratic-quadratic 
combinations  are  used  to  approximate  the  distributions  of  displacements 
and  stresses  in  the  mixed  model.  All  of  the  computations  were  carried 
out  on  an  IBM  3081  in  double  precision  arithmetic. 
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The  cantilever  beam  shown  in  Figure  1  was  analyzed  using  5  8-node 
quadratic-quadratic  elements  using  the  assumptions  of  plane  stress.  The 
material  properties  and  geometry  used  are  given  in  the  figure.  The 
analysis  was  completed  using  16  equal  load  steps  and  a  tolerance  of  E  = 
0.001.  The  displacements  along  the  line  x  =  0.0  were  specified  to  be 
zero. 

Figure  2  contains  the  plot  of  nondimensionalized  tip  displacement 
vs.  applied  load  for  the  linear  analysis  and  nonlinear  analysis  as 
computed  by  the  mixed  and  displacement  finite  element  models.  The  mixed 
model  gives  displacements  that  are  larger  than  those  of  the  displacement 
model,  with  a  maximum  difference  of  2.8  percent  measured  at  the  final 
load  level.  Although  the  results  of  the  displacement  model  are  in 
excellent  agreement  with  the  analytical  solution  reported  by  Holden  • 
[24],  it  should  be  mentioned  that  Holden's  solution  uses  the  Euler- 
Bernoulli  beam  theory  (i.e.,  does  not  even  account  for  the  transverse 
shear  strain).  However,  for  a  beam  of  the  dimensions  used  in  this 
example,  the  shear  deformation  effect  is  undoubtedly  quite  small,  and  it 
remains  that  the  mixed  element  model  yields  slightly  larger 
displacements  than  does  the  displacement  model  for  this  example. 

It  is  also  of  interest  to  compare  the  Cauchy  stresses  determined 
from  the  two  formulations.  In  the  displacement  model,  it  is  necessary 
to  compute  the  Almansi  strain  tensor,  and  then  use  the  constitutive 
relations  of  the  material,  as  described  in  Eq.  (26).  In  the  mixed 
model,  we  may  simply  interpolate  the  values  of  the  nodal  stresses  using 


Thickness,  t  =  1  in. 


Figure  1  Geometry,  material  properties,  and  finite  element 
mesh  used  for  uniformly  loaded  cantilever  beam 
example. 


Load  parameter,  K 
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Figure  2  Tip  displacement  of  uniformly  loaded  cantilever  (K  =  PL  /El) 


Distance  along  the  beam,  x  (in.) 


Figure  3  Bending  stress  along  upper  Gauss  points  of  cantilever  beam 


the  original  form  of  the  stress  approximation,  i.e. 

Tij(xi.  x^)  = 

To  compare  the  Cauchy  stresses  using  the  two  formulations,  we 
computed  the  axial  stress  at  the  location  of  the  2x2  Gauss  points 
along  the  upper  half  of  the  beam.  The  displacement  components  from  the 
two  models  are  not  exactly  the  same  (the  actual  positions  of  the  Gauss 
points  vary  somewhat),  but  this  difference  is  very  small.  The 
comparison  of  the  axial  stress  components  is  shown  in  Figure  3  for  the 
two  different  models  at  the  final  load  level.  The  agreement  is 
excellent,  with  the  mixed  results  giving  a  maximum  stress  that  is  4.8 
percent  higher  than  the  maximum  stress  computed  using  the  displacement 
model.  All  other  computed  values  are  in  much  closer  agreement  for  the 
two  models. 

Although  the  displacements  in  the  y-direction  are  in  very  good 
agreement  for  the  two  formulations,  the  difference  in  the  displacement 
gradients  can  be  quite  large  and  can  result  in  significant  errors  when 
computing  the  Cauchy  stresses  from  the  Almansi  strains  using  the 
displacements  from  the  mixed  formulation.  For  example,  considering  the 
lower  left  Gauss  point  of  the  rightmost  element,  we  have  the  following 
values  at  a  given  load  step: 

Displacement  Formulation  Mixed  Formulation 
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0.003 

0.007 

^y 

0.082 

0.082 

^x 

-0.082 

-0.082 

'',y 

0.003 

0.007 

^xx 

-1.95 

51.20 

The  stress  computed  from  the  displacement  model  is  in  very  good 
agreement  with  the  stress  computed  from  ’■he  mixed  model  using  nodal 


stress  interpolation,  but  the  stress  computed  from  the  Almansi  strains 
using  the  displacements  from  the  mixed  model  is  very  much  in  error.  It 
appears  from  the  results  of  this  example  that  the  nodal  displacements  of 
the  mixed  model  should  not  be  used  at  any  point  of  the  analysis  to 
compute  the  Cauchy  stress  components,  and  the  nodal  stress  values  should 
be  used  instead. 


Rigid  Pin  Examples 

The  contact  algorithm  proposed  by  Wilkenson  and  Rahman  was 
implemented  using  the  three  basic  iterations  of  load,  contact,  and 
friction  using  a  geometrically  nonlinear  formulation  along  with  mixed 
finite  elements.  Here  the  results  of  several  example  problems  involving 
contact  between  an  elastic  body  and  a  rigid  pin  are  presented  not  only 
to  demonstrate  the  accuracy  of  the  algorithm  but  also  to  highlight  the 
effectiveness  of  having  stress  as  a  nodal  variable.  The  results  of  the 
example  problems  are  compared  with  available  analytical  and  numerical 
solutions. 

In  the  first  example  we  consider  an  infinitely  long  cylinder  of 
radius  r  =  1  resting  on  a  rigid  plane  and  under  a  uniform  line  load. 

This  problem  was  modeled  using  the  assumptions  of  plane  strain.  One 
quarter  of  the  circular  domain  was  used  to  model  the  problem,  and  was 
approximated  by  84  linear-linear  elements.  The  load  was  assumed  to  act 
at  the  center  of  the  cylinder  and  was  applied  in  12  increments,  with  the 
initial  increments  smaller  than  the  later  increments.  The  problem  was 
modeled  by  assuming  that  the  cylinder  was  in  contact  with  a  rigid  pin  of 
very  large  radius  (R  =  1000)  to  model  the  rigid  plane.  A  tolerance  of  E 
=  0.001  was  used  for  the  equilibrium  iterations.  The  modulus  of 
elasticity  used  was  21,000  ksi  and  the  Poisson's  ratio  used  was  0.3. 


The  results  of  the  analysis  are  shown  in  Figures  4  and  5.  The 
numerical  results  are  compared  with  the  Hertz  analytical  solution. 

Figure  4  shows  the  contact  pressure  distribution  plotted  against  the 
distance  from  the  original  point  of  contact.  Figure  5  shows  the  load 
plotted  against  the  total  contact  area,  the  data  points  of  which  can 
only  be  determined  when  each  successive  node  comes  in  contact  with  the 
pin.  The  results  from  the  contact  algorithm  appear  to  be  quite  good. 

The  second  example  considers  a  thin,  orthotropic,  pin  loaded  plate 
with  a  hole  of  radius  r  under  a  uniform,  in-plane  load  as  shown  in 
Figure  6.  Due  to  symmetry,  one  half  of  the  plate  was  modeled  using  124 
linear-linear  elements  with  156  nodes  for  a  total  of  780  degrees  of 
freedom.  The  material  properties  given  in  the  figure  are  averaged 
properties  from  a  number  of  species  of  wood.  The  pin  was  assumed  to  be 
rigid  and  of  radius  R.  The  plate  was  loaded  to  a  final  load  of  400 
pounds  per  inch  of  plate  thickness,  and  was  applied  in  18  increments.  A 
coefficient  of  friction  of  0.7  was  assumed  at  all  points  of  contact 
between  the  plate  and  the  pin  for  all  load  levels,  and  is  a  typical 
value  for  wood  on  steel.  No  equilibrium  iterations  were  performed  for 
this  problem. 

Figure  7  shows  the  radial  stress  distribution  as  a  function  of  the 
angular  position  around  the  pin  for  the  nodes  that  have  come  in  contact 
at  the  final  load  step.  These  results  are  compared  with  the  results 
obtained  by  Wilkenson  [181  using  a  finer  mesh  (385  nodes)  and  quadratic 
displacement  elements.  The  comparison  is  very  good. 
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